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ABSTRACT 

The usefulness of angle-action variables in galaxy dynamics is well known, but their 
use is limited due to the difficulty of their calculation in realistic galaxy potentials. 
Here we present a method for estimating angle-action variables in a realistic Milky 
Way axisymmetric potential by locally fitting a Stackel potential over the region an 
orbit probes. The quality of the method is assessed by comparison with other known 
methods for estimating angle-action variables of a range of disc and halo-type orbits. 
We conclude by projecting the Geneva- Copenhagen survey into angle-action space. 

Key words: methods: numerical - The Galaxy: kinematics and dynamics - galaxies: 
kinematics and dynamics - The Galaxy: solar neighbourhood - The Galaxy: structure 
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1 INTRODUCTION 

In the study of dynamical systems it is becoming increas- 
ingly important to be able to process and understand large 
multi-dimensional data sets efficiently. The stars in our own 
galaxy, the Milky Way, are being increasingly observed in 
full six-dimensional phase-space through the combination 
of astrometry and radial velocity measurements. Full 6D 
phase-space information is currently available for stars in 
the solar neighbo urhood from the Geneva- Copenhagen and 
RAVE surveys (jNordstrom et al.ll2004l : IZwitter et al]|2008l ) 
and this is to be greatly expanded o n by the future space 
mission Gaia ([Perrvman et al.ll200lh . Beyond our Galaxy, 
the advent of integr al-field spectroscop y has led to projects 
such as SAURON (jBacon et al.ll200lL and subsequent pa- 
pers) , which mapped the kinematics of a representative sam- 
ple of 72 nearby ell iptical and spiral ga laxies, and subse- 
quently ATLAS3D ([Cappellari et al" l l201lL and subsequent 
papers), which combined SAURON observations with CO 
and HI observations to study the kinematics of a complete 
volume-limited sample of 260 local early-type galaxies. Ob- 
servational data is often understood by performing large N- 
body simulations. Whilst such models are straightforward 
to produce, the configurations of the models are difficult to 
control and characterise. Schwarzschild modelling offers an 
improvement on this by describing the configuration of a 
model by a weighted set of orbits. However, this approach 
is not the most natural as each orbit is characterised by 
its initial phase-space coordinates. It is necessary that tech- 
niques are developed which can simplify both observational 
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and simulation data without losing the richness of the phase- 
space information. 

Angle-action variables are a set of canonical coordinates 
which can be used to express the equations of motion in a 
trivial form: the actions are integrals of the motion whilst 
the angles increase linearly with time. Such a formulation 
instantly reduces the complexity of any dynamical data set 
by reducing the six phase-space dimensions to three angle 
coordinates. Angle-action variables can be defined for any 
quasi-periodic orbit. Initially introduced to study celestial 
mechanics, angle-action variables now have great potential 
for galaxy dynamics due to their attractive properties. For 
instance, the Jeans theorem states that the arguments for 
the distribution function of a steady-state galaxy must be in- 
tegrals of the motion, and it is particularly convenient to use 
as integrals the actions as (i) they are adiabatic invariants, 
(ii) the zero-point of an action is well defined and (iii) the 
range of values an action may take is independent of the 
other actions. The angle-action variables also provide a basis 
for the development of a perturbative solution to the equa- 
tions of motion (see iBinnev fc Tremaind 120081 for a much 
fuller discussion of the merits of angle-action variables). 

The increasing evidence of s ubstructure within the stel- 
lar halo of the Milky Way (e.g. iBelokurov et al.ll20Qq ) has 
led many authors to consider the use of angle-action vari- 
ables when hunting for and understanding the formation 
of structure within phase -mixed data sets. For example, 
iMcMillan fe Binnevl (2008) studied a simulation of a self- 
gravitating satellite in a realistic Galaxy potential in angle- 
action space. Though the stars became well phase-mixed, 
the action space still showed considerable structure and 
through the use of angle- variable diagnostics the Galaxy po- 
tential and history of the satellite could be reconstructed. 
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Similarly, ISellwoodl (|2010h and lMcMilknl (|2011aT ) have used 
the Geneva- Copenhagen survey to analyse stars in the solar 
neighbourhood in angle space and showed that the Hyades 
moving group may be due to a recent inner or outer Lind- 
blad resonance. The study of tidal streams has a natural 
expression in angle-action variables. Under certain condi- 
tions the dimensionality of the stream may be reduced to 
one a s the stream st r etches out in a single angle coordi- 
nate (|Tremaindll999h . lEvre fc Binnevl (|201lh showed that 
the path of a stream can be reconstructed far more reli- 
ably in angle-action space than by incorrectly assuming that 
streams delineate the orbits of their progenitors. 

Despite the aforementioned advantages, angle-action 
variables remain awkward to work with in practical appli- 
cations due to the difficulty of their calculation in a gen- 
eral potential. They are easily calculated when the potential 
is spherical and with more work can be analytically calcu- 
lated when the potential is of Stackel form, but neither of 
these approaches is satisfactory when working with realis- 
tic galaxy potentials, because such potentials do not satisfy 
these conditions. The development of methods to estimate 
angle-action variables in a general potential is crucial if we 
are to benefit from the advantages of angle-action variables 
and the wealth of techniques which utilise them. 

In this paper we present a method for estimating angle- 
action variables in a general axisymmetric potential. The 
method proceeds by fitting a Stackel potential locally to the 
region of the potential a given orbit explores, thus enabling 
us to calculate analytically the actions and angles in this fit- 
ted Stackel potential. In Section [2] we give a brief overview 
of the determination of angle-action variables in an axisym- 
metric Stackel potential and then in Section [3] present the 
method for locally fitting such a potential to any axisym- 
metric potential. The results of the method are examined 
by analysing artificial data in Section [5] and then these re- 
sults are compared to other methods in Section [6] Finally 
we demonstrate the practical application of the method by 
inspecting the Geneva- Copenhagen Survey in angle-action 
space in Section [71 



2 ACTIONS AND ANGLES IN A STACKEL 
POTENTIAL 

The most general class of potentials in which we are able 
to calculate the angle-action variables analytically is that 
of Stackel potentials. In a confocal ellipsoidal coordinate 
system these potentials produce separable Hamilton- Jacobi 
equa tions. A full discussion of Stackel potentials is given 
m Ide Zeeuwl (1985). Here we limit the discussion to oblate 
axisymmetric Stackel potentials which are associated with 
prolate spheroidal coordinates (A, 0, v). A specific prolate 
spheroidal coordinate system is defined by two constants 
(a,c). These coordinates are related to cylindrical polar co- 
ordinates (R, <fi, z) by 



R 2 



+ 



1, 
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where A and v are the roots of r such that c 2 ^ v ^ a 2 ^ A. 
Surfaces of constant A are prolate spheroids and surfaces of 
constant v are two-sheeted hyperboloids of revolution which 
intersect the spheroids orthogonally. A potential, is of 



Stackel form in a particular prolate spheroidal coordinate 
system if 
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&s is fully defined by a single function f(r). A single func- 
tion may be used as A and v take different ranges of values 
except at A = a 2 , v — a 2 , where we require / to be contin- 
uous so the potential remains finite. As in all axisymmetric 
potentials, the energy, E, and ^-component of the angular 
momentum, L z are isolating integrals. In a Stackel poten- 
tial we are in the fortunate position of being able to find 
analytically a third isolating integral, 
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Therefore, given a Cartesian phase-space point (cc,v), we 
can find the three isolating integrals, I = (E,L z ,l3), us- 
ing the coordinate transformation and calculating A from 
v. Using the three isolating integrals we can write the first 
integrals of motion as 

2(T _ aK= £_______ + __ (4) 

where p T is the momentum conjugate to r = A, v. We define 
the action variables, J\ and J u , as 



J T = — (p p T dr, 



(5) 



where the integration is over all values of r for which p 2 ^ 0. 
As p T = Pt(t, E, 12,13) the actions are solely functions of the 
isolating integrals and thus constants of the motion. The 
third action, J^, is simply L z . The actions give an absolute 
measure of the extent of the oscillations of the orbit in each 
of the coordinates. At large radii the prolate spheroidal co- 
ordinate system becomes spherical such that A « R 2 + z 2 . 
Therefore we can think of J\ as a measure of the radial os- 
cillations. The v coordinate increases as we move away from 
the z — plane so we may think of J v as a measure of the 
vertical oscillations. 

The corresponding angle coordinates, T and 0$, are 
calculated by the introduction of the generating function, 
S(X, <fi, v, J\, L z , J u ) for the canonical transformation from 
(\,(\),v,p\,p v ,L z ) to {0\,0^,0 v ,J\,L z ,J v ). The angles are 
found by differentiating the generating function with respect 
to the respective action such that 
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for t — A, z/; 
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(6) 



A full list of formulae, as well as a discussion of how to per- 
form the quadratures numerically, is given in Appendix [A] 



3 FITTING A POTENTIAL WITH A 
STACKEL POTENTIAL 

Given the ease with which we can calculate actions and 
angles in a Stackel potential, it seems sensible to investi- 
gate how well a Stackel potential can fit a Galaxy model so 
that we may calculate the actions and angles in this best-fit 
potential. It has been known for some time that Stackel 
potentials do not give a good fit to the potential of the 
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Galaxy globally due to the rigid conditions they must fulfil. 
IDeionghe &; de Zeeuwl (| 19881 ) outline a method for fitting 
an axisymmetric potential with a Stackel potential, which 
can be applied both globally and locally. These authors pro- 
duced global fits for the B ahcall-Schmidt-Soneira Galaxy 
model dBahcall et al. 1982h with errors nowhere exceeding 
3 per cent, and I Jaseviciusl jl994) carried out a similar anal- 
ysis on a broader range of Milky Way potential models with 
similar results. As expe cted, the fits are worst in the central 
0.5 kpc of the Galaxy. lDe~Bruvne et all tOQO) sought to fit 



axisymmetric potentials locally using a set of Stackel po- 
tentials in order to calculate the third isolating integral, I3. 
When applied to a Miyamoto- Nagai potential, ^3 was found 
to vary by approximately 10 per cen t along an orbit. Here 
we fo llow the method presented by IDeionghe &; de Zeeuwl 
(| 19881 ). 

Suppose we have an axisymmetric potential z) 
that we wish to fit by a Stackel potential <J>fit- We begin 
by choosing a prolate spheroidal coordinate system by spec- 
ifying (a,c). The coordinate system is fully specified by the 



combination 



") so we are free to set c =1, which re- 



duces numerical difficulties. We determine a by using a prop - 
erty of an axisymmetric Stackel potential (|de Zeeuwl 1 19841 ). 
It follows from equation ([2]) that for a Stackel potential, $5, 



-[(\-v)<S>s] =0. 
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Therefore, for a general potential <I> we use this equation as 
a definition for the coordinate system. Using expressions for 
R and z as a function of (A, v), 

2 _ (A-a 2 )(^-a 2 ) 
u ~ 2 _ 2 ' 

_ (A-c^-c 2 ) (8) 
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we find this gives an estimate for a at a point (R, z 

dR 2 dz- 
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We calculate a sufficiently accurate value of a by evalu- 
ating this expression at multiple positions along the orbit 
and averaging. With this choice of a we transform &(R, z) 
to <&(A, v) and specify the fitting region: A_ ^ A ^ A+, 
v- ^ v ^ z/+ . A global fit corresponds to v- — c 2 , v+ — 
A- = a 2 , A+ = 00. We also define the auxiliary function 



x (A,i/) = -(A-i/)*(A,i/). 



(10) 



If the potential <I> is of Stackel form, this auxiliary function is 
simply x(A, v) — /(A) — f(y). We seek the function / which 
makes <3>fit most like <I> by minimising the square difference 
of the potential auxiliary function and the fitting potential 
auxiliary function, %fi t , over the fit region. Therefore, we 
minimise the functional 



F[f] = dX dv 

J\_ Jv_ 



(11) 

where A (A) and N(y) are weighting functions allowing us to 
acquire a better fit in certain areas. These functions must 
be finite when integrated over the fitting region. We choose 
the normalised weighting functions 

a(a) = 4A- 5 (a: 4 -a; 4 )-\ n(u) = {v+ - V -)~\ (12) 



This choice of weighting functions gives preferential weight 
to smaller values of A where the potential is harder to fit. 
Analytic minimisation of the functional F results in a best 
fit function 



/(A) = x(A) - \x, fiy) = -xiy) + \x (13) 



where 



X(A) = /" + dv X {\v)N{v) 
XW) = j + dA X (A,^)A(A) (14) 
X= f X+ ' dXdu X {X,u)\{X)N{u). 

Jv_ J\_ 

The derivation of these equations is given in Appendix [B] 
The quality of the fit we have achieved is then measured by 
F[f]- 

4 PROCEDURE 

Combining the above two sections we can estimate the ac- 
tions and angles of a phase point (x,v) by first fitting a 
Stackel potential to the given potential over the region the 
orbit probes and then calculating the angle-action variables 
in this fitted potential. Therefore, given a point (x,v) we 
follow this procedure: 

(i) We begin by calculating the z-component of the angular 
momentum, L z , and the energy, E, in the 'true' potential. 

(ii) We then integrate the orbit in the 'true' potential. We 
use the initial time-steps of the orbit integration to find the 
best-fit coordinate system: at several points along the orbit 
we evaluate equation ((9]) and average to find a sufficiently 
accurate value for a. With the coordinate system found, we 
continue integrating to find the edges of the orbit A+,A_ 
and zy + , which define the fitting region. The edges of the 
orbit are given approximately by the points where f = 
in the best-fit ellipsoidal coordinate system. The minimum 
and maximum r edges are distinguished by inspecting the 
sign of t. For all realistic potentials every orbit crosses the 
z = plane so we set v — c 2 . 

(iii) We can now find a best-fit Stackel potential over this 
region. Using equation ([T3|) we tabulate /(A) and f(y) for 
40 points in (A_,A+) and (z/_,z/+) respectively so that we 
may interpolate these smooth functions. Any call outside 
the ranges is calculated fully using equation ([T3|) with a full 
re-computation of x( T )- 

(iv) With the best fit potential now calculated, we find I3 
using equation (|3]) for three points on the boundary of the 
orbit (on the minimum A edge, the maximum A edge and 
the maximum v edge) and take an average. We have already 
found these three points when determining the edges of the 
orbit so this choice involves minimum additional computa- 
tional effort and provides a fair estimate for ^3 over a large 
region of the orbit. However, this choice of ^3 can lead to 
the initial phase-space point (x,v) being forbidden. There- 
fore, with this choice of ^3 we check whether p 2 (y) > 
and p 2 (X) > for the initial phase-space point using equa- 
tion and if not then we calculate ^3 from equation Q 
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Figure 1. Contours of ln($(i?, z)/$(0, 0)) for McMillan's best-fit 
Milky Way potential. The contours are increasing from the centre 
in equally space units of 0.15 with the central contour at —0.15. 



using only the initial phase-space point. This procedure re- 
duces the numerical noise around the turning points, partic- 
ularly in R. 

(v) With the three isolating integrals calculated, we are in a 
position to estimate the actions and angles using the method 
outlined in Section [2] The limits of the orbit are redeter- 
mined by finding from equation (4J the points where p 2 T — 
and are not given by r±. 

Table [T] quantifies the efficiency of this procedure. 



using the 'torus machine ' (|McMillan fc Binnevl hoOS). Or- 
bital tori are three-dimensional surfaces characterised by the 
three actions J = (Jr, L z , J z ) obtained as the images of an- 
alytic tori under a canonical transformation. The strength 
of the torus machine lies in constructing a torus given a set 
of actions, J, such that the phase-space coordinates, (sc, v), 
may be obtained as functions of the angles, 0, over the sur- 
face of the torus. Therefore, a simple test for the Stackel po- 
tential fitting procedure is to produce a list of phase space 
coordinates with fixed actions but randomly chosen angles 
using the torus machine. The success of the method is then 
measured by how accurately the angle-action variables can 
be reproduced. We note that the canonical transformation 
produced by the torus machine maps Jr into J\ and J z into 
J v . From now on we will use the more intuitive notation for 
the actions, Jr and J z , and similarly for the angles, Or and 

e z . 

The errors in the actions of a given torus from the torus 
machine may be estimated from the residuals of the Hamil- 
tonian over the surface of the torus. The error in the Hamil- 
tonian, AH, is related directly to the error in one of the 
actions by 



AH 



DH 
dJ 



A J — Q A J 



(15) 



where we find the frequency Q directly from the torus ma- 
chine. Assuming the errors in Jr and J z are approximately 
equal and uncorrelated, the error in the actions may be es- 
timated as 



A J : 



AH 



(16) 



The true angles of an orbit in a potential increase linearly 
with time. The errors in the torus angles are estimated by 
the residuals of the angles away from this expected straight 
line. Clearly we require these errors to be smaller than the 
errors from the Stackel fitting procedure in order to state 
anything meaningful about the systematic errors from our 
method. 



5 APPLICATION 

We now investigate how successful the above routine is in 
calculating the action-angle variables in a general axisym- 
metric potential. To demonstrate the applicability of the 
met hod to data we ch oose a realistic Milky Way potential 
from McMillan (2011a). This potential consists of two expo- 
nential discs for the thick and thin discs of the Galaxy and 
two spheroids for the bulge and dark matter halo. We select 
the 'best' model from this paper. The equipotential contours 
for this model are plotted in Figure [T] It is clear that as we 
move out from the centre, the contours become more circu- 
lar so we anticipate that they are better fit by surfaces of 
constant A and v. Therefore, we expect more accurate esti- 
mates of the angle-action variables for orbits at larger radii. 
Also orbits that probe a large range of R and/or z should 
have less accurate angle-action variable estimates as these 
orbits probe a large range of curvature of the equipotential 
contours. Therefore, we expect the method to work best for 
small J\ and J u but large L z . 

We assess the validity of the method by comparing 
the results with the 'exact' angle- action variables calculated 



5.1 Single Torus 

Here we discuss the results of applying the procedure to 
10000 randomly generated points from the toru^l J — 
(J R , L z , J z ) = (0.078, 1.9, 0.097) kpc 2 Myr _1 . This torus was 
chosen to be representative of the actions of a disc star 
in the solar neighbourhood. For this torus the errors in 
the actions and angles are A J/ J = 0.01 per cent and 
(A0 R , A6> , A0 Z ) = (1.0,0.2,1.0) x 10" 5 rad. The orbit in 
the (R,z) plane is shown in Fig. [2] This orbit has apses at 
R ~ (6.5, 10.5) kpc and z ma x ~ 2.8 kpc. Also shown in the 
figure are the curves defining the fit region and equipotential 
contours for McMillan's best-fit potential. 

The residuals in the fitted potential over the fitting re- 
gion defined in Fig. [2] are plotted in Fig. [3] Everywhere 
within the fitting region the error in the potential is less 
than 0.2 per cent of the maximum difference in the poten- 
tial across the fitting region. We note here that a good fit for 



1 Throughout this paper the actions are stated in units of 
kpc 2 Myr _1 = 977.8 kpc km s -1 




Angles and actions in an axisymmetric potential 5 



Figure 2. Fit region for a single orbit - the blue 
line shows the orbit with actions J = (Jr,L z ,J z ) = 
(0.078, 1.9, 0.097) kpc 2 Myr- 1 . The black lines are the equipoten- 
tial contours of McMillan's best-fit Milky Way potential. The red 
lines show the lines of constant A and u which define the region 
over which the potential is fitted. 



the potential does not necessarily correlate with an accurate 
calculation of the actions. Small changes in the potential can 
cause large changes in the motion of a particle so, whilst a 
good fit for the potential is necessary, we don't expect the 
errors in the actions to be of similar order. 

The 10000 phase-space points are shown in scatter plots 
of (R,9r) and (z,0 z ) in Fig. [H We can see that R and z 
are periodic in the angles. We define the zero-point of Or 
such that the radial periapsis and apoapsis correspond to 
Or = and Or = n respectively. Z is defined such that 
z = corresponds to Z = 0, tt and z = =bz ma x corresponds 
to Z = 7r/2, 37t/2. The zero-point of 0$ is defined such that 
0^ = 4> at periapsis. The spread of the z coordinates of the 
points at a given angle is much larger than the spread in the 
R coordinates. 

When the Stackel fitting method is applied to this 
set of phase-space points, we find that the root-mean- 
square (RMS) deviations of the actions, AJ, are given by 
AJr/Jr psj 4.9 per cent and AJ Z /J Z « 4.2 per cent. We 
also find that there is a very tight anticorrelation between 
Jr and J z . All phase-space points have the same energy as 
we are using the potential that was used to integrate the 
orbit to find the energy. Therefore, all the points lie along 
the intersection of the surface of constant energy with the 
(Jr, Jz) plane. If we overestimate Jr then we must under- 
estimate J z in order to have the correct energy. The RMS 
deviations in the angles for the 10000 phase-space points are 
(A0 R , A0+, A0 Z ) = (4.1, 1.1, 5.1) x 10~ 2 rad. 



5.1.1 Errors as a function of angle 

It is informative to investigate how the errors in the derived 
actions and angles vary with true angle around the torus. 
The derived actions are approximately independent of the 
true angles as they depend only on the path of the orbit, 
which is determined by the fitted potential and not the ini- 
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Figure 3. Filled contour plot of the percentage difference be- 
tween the best-fit Stackel potential and McMillan's best-fit Milky 
Way potential. ^ m i n and $ max give the values of the potential 
on the minimum and maximum A edges respectively. Also plot- 
ted in black are the curves of constant A and v which define 
the region over which the potential is fitted. This is the fit re- 
gion corresponding to the orbit shown in Fig. [2] with actions 
J = (J R ,L Z ,J Z ) = (0.078, 1.9, 0.097) kpc 2 Myr- 1 . We see that, 
within the fitting region, the difference between the fitted po- 
tential and the potential we are attempting to fit is less than 
0.2 per cent of the maximum potential difference across the fit- 
ting region. 



tial point on the orbit. Any small variation is due to the 
choice of ellipsoidal coordinate system and variations in the 
fitted potential. However, we find that the error in the de- 
rived angle varies with true angle. In Fig.[5]we plot the RMS 
errors in the angles binned as a function of true angle for 
both Or and Z . Maximum errors occur at the turning points 
in the (R, Or) and (z,6 z ) plots shown in Fig. [4] For the ra- 
dial and vertical angle the largest error occurs at apoapsis. 
In a Stackel potential the momenta, p T , depend on r and 
the isolating integrals, which once determined are taken to 
be constant. Therefore at a given location in the orbit the 
angle is solely a function of the position coordinates and the 
velocity information is essentially ignored. Around turning 
points in the orbit the velocity coordinates contain the ma- 
jority of the information whilst the position coordinates are 
changing very slowly. Therefore at turning points the errors 
in the angles are large as the angle coordinates are estimated 
using this reduced phase-space information. In general the 
errors in Z are larger than the errors in Or. 



5.2 Multiple Tori 

We have seen that the method gives reasonable esti- 
mates for the actions for a particular torus, but in or- 
der to use the method with confidence we need to see 
how the errors depend on the torus. Here we repeat 
the above procedure for a range of different tori which 
probe the different regions of the potential. We work 




Figure 4. Scatter plot of R against Or and z against Z for 10000 
randomly selected phase-space points from the torus detailed in 
Section 15.11 

with two groups of tori: those with low actions and 
torus machine errors less than A J/ J — 0.01 per cent and 
(A0 R , A6> , A6 z ) = (27.0,5.1,990) x 10" 6 rad and those 
with high actions and torus machine errors less than 
AJ/J = 1 per cent and (A0 R , A6> , A6 z ) = (2.0, 1.7, 1.7) x 
10 -2 rad. The low-action group consists of 100 tori 
with actions J R = (0.001, 0.005, 0.01, 0.05, 0.1) kpc 2 Myr -1 , 
J z = (0.001, 0.005, 0.01, 0.05, 0.1) kpc 2 Myr- 1 and L z = 
(1.0,2.0,3.0,4.0)kpc 2 Myr _1 . These tori probe the re- 
gion 3kpc < R < 22kpc, z < 5kpc and are cho- 
sen to be representative of disc-type tori. The high- 
action group consists of 36 tori with actions Jr = 
(0.5,1.0,5.0) kpc 2 Myr _1 , J z = (0.5, 1.0, 5.0) kp^Myr" 1 and 
L z = (1.0, 2.0, 3.0, 4.0) kpc 2 Myr _1 . These tori probe the re- 
gion 2kpc < R < 120 kpc, z < 100 kpc. We include the 
second group to demonstrate that the method can deal with 
orbits which deviate very far from the plane and probe a 
very large region of the potential. We would like to be able 
to apply the method to halo stars and tidal streams so it 
is important to understand the errors for these high-action 
tori. 



5.2.1 Actions 

As mentioned previously, we expect the errors in Jr and J z 
will be large when Jr/\L z \ and/or J Z /\L Z \ are large. In this 
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Figure 5. RMS error in the angles binned as a function of angle. 
The dashed line shows the total RMS error from all the points on 
the torus. For the bottom panel we have taken advantage of the 
symmetry in the z = plane and mapped the Z = (0, 2tt) interval 
onto 6 Z = (0, 7r) such that 6 Z = tt/2 corresponds to ±z max etc. 
The largest error occurs at the apses for both cases. 



regime the orbit probes a large central region of the potential 
so we anticipate the potential fit will be poorer. In Fig. [6] 
the RMS deviations in the actions for the complete orbit 
sample are plotted against the combination of the actions 
(Jr-\-Jz)/\L z \. We can see that, as anticipated, the absolute 
errors correlate with this action combination. In fact, the 
correlation is much tighter than the individual correlations 
with Jr/\L z \ and J Z /\L Z \, so the errors in the method are 
dependent on the sum of the actions (Jr + J z ). It is this 
measure which tells us how much an orbit strays from a 
circular orbit and thus how much of the potential it explores. 

We also note from Fig. [6] that at a given value of (Jr + 
J Z )/\L Z \ the errors in Jr and J z are of similar magnitudes. 
As explained above, the errors in Jr and J z compensate for 
each other to recover the correct energy. In Fig [71 we plot 
this correlation between the RMS errors in Jr and J z . A 
consequence of this tight correlation is that when one action 
is much greater than the other, the relative error in the 
smaller action will be much greater than the relative error 
in the larger action. However, it is worth noting that the 
absolute error is far more important than the relative error. 
Given a distribution function for a steady-state galaxy, /(J), 
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the absolute error in / is given by 
(A/) 2 



X^df df 



(17) 



where cov(X, Y) is the covariance between variables X and 
Y. In the case of uncorrelated errors between the actions 
this simply becomes 



(A/) 2 



\dJi 



(18) 



The distribution function for the Milky Way is approxi- 
mately exponential in the actions (|Binnevll2010l ): 



(19) 



where a% is independent of the action Ji. Therefore the ab- 
solute error in / is given by 



(a/) 2 = Y^u^Ji. 



(20) 



Similarly the relative error in the distribution function is 
given by 



£( 



d In/ AJ t 
d In Ji Ji 



) 2 = Y / (a i AJ l f. 



(21) 



Both the absolute and relative error in the distribution func- 
tion are determined by the absolute errors in the actions, so 
we need not be overly concerned that the relative error in 
one action is much larger than the relative error in another. 

From the relationship illustrated in Fig. [6] we can esti- 
mate the error in a given estimate of Jr and J z . Performing 
a linear fit to both sets of data points independently, we find 
that for (J R + J Z )/\L Z \ < 10 a good fit for the RMS errors 
in both Jr and J z is given by 
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Figure 6. Absolute RMS deviations of the actions. The blue cir- 
cles are data for which the relative errors in the actions from the 
torus machine are less than 0.01 per cent and the red squares are 
those with errors less than 1 per cent. The errors correlate loosely 
with both Jr and J z separately, but there is a much tighter cor- 
relation between the errors and (Jr + J z ). 



(22) 



The errors in the actions have a weak dependence on L z as 
orbits with higher L z explore regions of the potential which 
are more spherical and hence easier to fit with a Stackel 
potential (see Fig. [I]). 



5.2.2 Angles 

We present the RMS deviations in the angles for the 100 tori 
in the low-action group in Fig. [8] In general the errors in the 
angles are larger than the errors in the actions. The calcu- 
lation of an action involves only a single integral whereas 
the corresponding angle calculation involves nine integrals. 
Each integral folds in more error from the fitting method 
so we expect the errors in the angle variables to be signif- 
icantly larger than the errors in the actions. From Fig. [8] 
we see that the errors in the angles correlate with the rel- 
ative error in the actions. The error in 0^ has been plotted 
against AJr/L z whilst the other two angles have been plot- 
ted against the relative error in their respective action. As 
the errors in Jr and J z are approximately equal (Fig. [7} 
the three sets of points are essentially 6i against AJ/Ji. As 
with the errors in the actions, we may estimate the error in 
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Figure 7. RMS deviations in Jr and J z . The blue circles are 
data for which the relative errors in the actions from the torus 
machine are less than 0.01 per cent and the red squares are those 
with errors less than 1 per cent. 
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Figure 8. RMS deviations in the angles for the 100 low-action 
tori as a function of the relative error in their respective action. 
For Off, we have plotted the error against AJr/L z . 



Figure 9. RMS deviations of J z for the Stackel fitting method 
(blue circles) and the RMS deviations in the spherical vertical 
action (L — \L Z \) (red squares) for the same sample of orbits. 



a given calculation of the angles by fitting the data in Fig. [8] 
We find the errors are approximately given by 

/AJ\0-75 /AJ\0-75 /AJ\0-5 

A^(-) ,A*.«( X ) ,A*,«(— ) . 

(23) 



6 COMPARISON WITH OTHER METHODS 

6.1 Total Angular Momentum 

Some authors have hunted for structure in the distribution 
of stars in spaces defined by phase space funct i ons ot her than 
the actions. For example, iHelmi &; de Zeeuwl (|2000h use the 
set of variables (E, L z , L), where L — \x x v\ is the total an- 
gular momentum, to attempt to find substructure within nu- 
merical simulation s of disrupted satellite galaxies, whereas 
Helmi et al. (2006) considered the 'APL space' of apocen- 
tre, pericentre and z-component of the angular momentum 
in order to identify signatures of past accretion events in 
the Geneva- Copen hagen Survey of the solar neighbourhood 
([Nordstrom et al.l |2004 ) . The total angular momentum is 
only conserved when we are considering spherical poten- 
tials. In a spherical potential the vertical action is simply 
J z = L — \L Z \. Here we investigate how much better we are 
doing when we estimate the vertical action using a Stackel fit 
than if we simply use L. Fig.[9]shows the absolute RMS error 
in the vertical action for the 100 low-action tori taken from 
the lower panel of Fig. [6] along with RMS error in the spher- 
ical vertical action, (L — \L Z \). The Stackel fitting method 
gives approximately two orders of magnitude improvement 
in the vertical action error compared to simply using L. 

6.2 Adiabatic Approximation 

The adiabatic approximation provides an alternative 
method for e stimating the actions. In its simplest form 
(jBinnevI [2010) this approximation assumes that the mo- 
tion in the plane is unaffected by the motion perpen- 
dicular to the plane. The absence of energy transfer be- 



tween the radial and vertical motion leads to an underes- 
timate of the centrifugal potential for the radial motion and 
hence an underestimate of the maximum radius of an orbit. 
iBinnev &; McMillan! (|201lh attempted to resolve this issue 
by repl acing L z with (\L Z \ + J z ) i n the effective radial po- 
tential. Schonri ch &; Binne vl (|2012h improved on this by in- 
cluding a correction to the radial energy due to the changes 
in the vertical energy along an orbit. It is this final approach 

that we test here. 

Following Sch onrich fe Binnevl |2012) we assume that 
the vertical motion at a given radius, i?o, is governed by the 
potential ^ z (z) = <f>(Ro, z) - <&(Ro, 0) such that the vertical 
energy, E z , is 

E z = ±v 2 z +* x (z). (24) 

Then the vertical action is estimated to be 

J z — — [ ZmaX dzv z , (25) 
^ Jo 

where z ma x is the point where the vertical velocity, v z , is 
zero. By linear interpolation we may reverse this calculation 
such that, for a given pair of J z and i?o, we may calculate 
E Z (J Z , Rq). Over the course of an orbit we take J z to be 
constant but the vertical energy will be changing as the orbit 
explores different radii. For overall energy conservation this 
energy must be transferred from the vertical motion into the 
radial motion. Therefore, the radial motion is governed by 
the one-dimensional potential 

L 2 

*r(R) = &(R,0) + ^+E Z (J Z ,R) - E z (J z ,R g ), (26) 

where R g is the guiding-centre radius (the radius of a cir- 
cular orbit with z-component of angular momentum L z ). 
Using this potential we estimate the radial action as 

1 f Ra 

Jr = - dRv R ,. (27) 

J R p 

where R p and R a are the points where the radial velocity, 
vr, is zero. 

This method is clearly much simpler than the full 
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Stackel potential fitting method and for a single phase-space 
point is approximately two to twenty times faster. Here we 
will see whether the accuracy of the Stackel potential fitting 
method justifies its extra expenditure. In Fig.[Tn]we have re- 
plotted the 100 low-action tori data from Fig. [6] along with 
the equivalent calculation using the adiabatic approxima- 
tion. We find that for low total action, (Jr + J z ), the Stackel 
fitting method is about two orders of magnitude more ac- 
curate than the adiabatic approximation, whilst at the high 
total action end this improvement is reduced to less than 
one order of magnitude. It is also worth noting that the rel- 
ative errors for the adiabatic approximation can be as much 
as order one for orbits with large J z - the adiabatic approxi- 
mation performs best for orbits that stay near the plane. In 
Table [1] we compare the errors in the actions and the time 
taken to estimate the actions for 10000 phase-space points 
for the Stackel fitting method and the adiabatic approxi- 
mation. We see that for near-circular near-planar orbits we 
can achieve an order of two magnitude decrease in the er- 
rors of the actions using the Stackel fitting method whilst 
only incurring an additional cost of doubling the computa- 
tional time. However, for orbits with much greater actions 
we can only achieve a single order of magnitude decrease in 
the error for a fivefold increase in the computational time. 



6.3 Iterative torus machine 



iMcMillan fe Binnevl (2008) used the torus machine itera- 
tively to find the actions and angles of a given phase-space 
point. This typically involved around 20 torus fits per phase- 
space point and relied on a good initial guess for the actions 
and angles for fast convergence. Such a method is poten- 
tially the most accurate way to determine the angle-action 
variables but also the most costly. The authors report that 
it takes around 15s to perform the iterative procedure on a 
single phase-space point. The alternative method we present 
in this paper could complement this method as it provides 
an accurate initial guess for the actions and angles which 
may be fed into the torus machine. We are using the torus 
machine as a source of 'true' actions and angles in order 
to test the accuracy of the methods, so we are assuming 
that, given enough iterations, this approach will give near- 
perfect results. Therefore, we do not explore the results of 
this method. 



7 GENEVA-COPENHAGEN SURVEY 

Now that we have understood the systematic errors of the 
method we ap ply it to real data. The Geneva-Copenhagen 
Survey (GCS) ([Nordstrom et alJl2004h is a sample of 16682 
nearby F and G stars, and is perhaps the best data set 
with full 6D phase-space information. It provides us with 
a platform to motivate a discussion of errors involved in a 
practical calculation. The GCS has been analysed by many 
authors and s pecificall y looked at in angle-action space by 
ISellwoodl (l2010h andlMcMillanl (l2011bh . From the table pro- 
duced bv iHolmberg et al.1 (|2009h we select the 13518 ob- 
jects for which we have the full 6D phase-space information. 
We correct the data for the solar velocity with respect to 
the local standard of rest as calculated bv ISchonrich et all 
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Figure 10. Absolute RMS deviations of the actions for the 
Stackel fitting method (blue circles) and the adiabatic approx- 
imation (red squares). We can see that in general the Stackel fit- 
ting method is an improvement on the adiabatic approximation. 
At the low total action ( Jr + J z ) end the improvement is around 
two orders of magnitude, whilst at the high total action end the 
improvement is reduced to less than one order of magnitude. 



i.e. (U,V,W)® = (11.1, 12.24, 7.2 5) kms" 



Using 



<&> . 

the 'best' model Galactic potential from McMillan] (|2011a ) 
sets the solar radius as R® =8.29 kpc and the velocity of the 
local standard of rest as vlsr = 239.1 kms -1 . The results 
of estimating the actions and angles using the Stackel fit- 
ting procedure are shown in Fig. 1111 These are v ery similar 
to th e equivalent plots from ISellwoodl {2010) and lMcMilhml 
(|2011bh . The shaded red region is inaccessible by stars which 
are at the solar position and have J z = 0. 

We see that the plot of Jr against L z has a markedly 
parabolic shape. The minimum of this parabola corresponds 
to the circular orbit at the solar radius. The edges of the 
parabola are defined by the limits of the survey. Stars with 
L z significantly different from the angular momentum of the 
local standard of rest, L z ^ = RqVlsr, require a sufficiently 
large radial action to bring them into the solar neighbour- 
hood. From the plot of 0$ against Or we see that the ma- 
jority of stars are at Or = 0, ty corresponding to the apses of 



2 We work in a right-handed Galactocentric Cartesian coordinate 
system with the positive x direction pointing towards the Galactic 
centre. 
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Table 1. Absolute RMS errors in the actions from the Stackel fitting method and the use of the adiabatic approximation along with the 
time taken to evaluate the actions of 10000 phase-space points for each method. 

Jr/\L z \ J z /\L z \ Aj| tack /AJAA Ajf tack /AJA A Timestack/s Time AA /s 



0.001 0.001 0.0087 0.12 10.1 4.5 

0.01 0.01 0.023 0.028 15.1 6.4 

0.1 0.1 0.16 0.16 44.6 8.3 
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Figure 11. Density contours for the actions and angles of 13518 
objects from the Geneva- Copenhagen Survey. The contours are 
linearly spaced. The red shaded region is inaccessible by stars 
which are at the solar position and have J z = 0. 



their radial motion. However, there is still a lot of structure 
in between these extrema: the peak at Or/tv psj 0.54 corre- 
sponds to the Hyades moving group. We do not plot the 
distribution in the vertical action and angle because these 
plots lack interesting features. 



7.1 Structure in the GCS 

Since iDehnenl (|l998l ) investigated the kinematics of the so- 
lar neighbourhood using data from the Hipparcos satellite, 
it has been known that when viewed in the (U, V) veloc- 
ity plane the local distribution of stars consists of a se- 
ries o££roup^_and clusters. These structures were classified 



by iFamaev et~ 



namical origin (|De 



(1200511 and are all thought to have a dy - 
■e Simone et al.ll2004l : lAntoia et aDl2Q10h . 



From the GCS sample we can identify several of the larger 



structures from the peaks in the (U, V) distribution. This is 
done by binning the stars in U and V and then perform- 
ing a wavelet transform. In Fig. [12] we show the results of 
the wavelet transform for structure on scales ~ 12kms _1 . 
The peaks in this plot are identified with the known groups, 
specifically the Hercules, H yades, Ple i ades, C oma Berenice 
and Sirius. As discussed by McMillan] (|2011bh for a star lo- 
cated at the solar position each point in the (U, V) plane rep- 
resents a unique point in the (Jr, L z , Or, 0$) space. Lines of 
constant L z and 0^ form an approximately Cartesian grid in 
the ([/, V) plane, whilst lines of constant Jr and Or form an 
approximately polar grid. Therefore, we can find a reduced 
set of angle-actions (excluding any vertical action and angle) 
for each of the peaks identified in Fig. [12] From inspecting 
the (U, V) plane we can see that the distribution deviates 
from axisymmetric equilibrium, so estimating the actions 
and angles assuming an axisymmetric potential is clearly an 
oversimplification. However, this is a necessary first step to 
create a basis on which we can include non- axisymmetric 
perturbations. 

We represent each group by a 6D phase-space point 
placed at the solar position with zero vertical velocity but 
U and V determined by the identification from Fig [12] and 
then estimate the actions and angles corresponding to this 
point, again using McMillan's 'best' potential. The results 
are shown in Table [2] This gives us an opportunity to discuss 
the various sources of errors in a realistic use of the method. 
There are two sources of error in this calculation - the 
systematic errors introduced by the Stacke l fitting method 
and the errors in the input coordinates. iHolmberg et al.l 
(2009|) estimated the space- velocity errors for each star as 
1. 5 km s -1 , which when c ombined with the errors estimated 
in lSchonrich et al. 1 (|2010h for the solar motion gives velocity 
errors of (A £7, AV) = (1.7, 1.6) kms -1 . The majority of this 
error arises from the uncertainty in the distances. We con- 
volve each peak position by these errors to calculate 10000 
Gaussianly distributed points in (£7, V) space about these 
peaks, and then estimate the errors in the output actions 
and angles by the RMS scatter in the resulting (J ,9) coor- 
dinates. These errors can then be compared and combined 
with the known systematic errors from Section [5.21 and are 
shown in Table [2] 

We find that for all the coordinates, the error in the data 
dominates t he systematic erro r introduced by the method. 
As noted by iMcMillanl (J2011b|) the relationship between er- 
rors in (x,v) and (J, 9) is non-trivial. Specifically at very 
low radial actions, any error in the velocity can introduce 
a 2tt error in the Or coordinate. We see this occurring for 
the Coma Berenice peak - the error in the Or coordinate is 
large as the peak is positioned very close to the origin of the 
(£7, V) plane. We also see that the errors in Jr and 0$ are 
of order one for Coma Berenice. 

We note that we have not included any error for the 
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Figure 12. Density contours for the (U, V) plane of 13518 ob- 
jects from the Geneva- Copenhagen Survey smoothed to a scale of 
~ 12 kms -1 via a stationary wavelet transform. The five major 
structures are identified by red points: 1. Hercules, 2. Hyades, 3. 
Pleiades, 4. Coma Berenice, 5. Sirius. 



size of the structures in phase-space, nor any error for the 
assumption that all the stars are situated at the solar posi- 
tion, nor any error in the choice of gravitational potential. 
Investigating the error in the actions due to the range of 
viable potentials for the Milky Way is beyond the scope of 
this paper. Even with this underestimated error, the errors 
from the data dominate the systematic errors. We conclude 
that, given the accuracy of the current data, the determina- 
tion of the angle-action coordinates using the Stackel fitting 
method is not limited by the well-understood systematic er- 
rors. 



7.2 Comparison with Mc Millan] (|2011bh 

iMcMillanl (|2011bl ) calculated the angles and actions of the 
GCS sample by using the torus machine iteratively. Here 
we compare the results to our own estimates of t he angle- 
actions. The potential used by McMillan (2011b) was the 
'convenient' potential detailed in McMillan! (|2011al ). which 
places the Sun at a galactocentric radius Rq = 8.5 kpc with 
a velocity of the local standard of rest vlsr = 244.5 kms -1 . 
We calculate the RMS deviations between McMillan's data 
and ours, and present the results in Table [3] We also show 
the expected RMS errors in the angles and actions estimated 
using the Stackel fitting method. The largest discrepancy 
between our data and McMillan's occurs for the Z coordi- 
nate. The expected error for this coordinate is also large, 
and very much larger than the actual error. This is because 
the sample is dominated by stars with low vertical action as 
the sample is situated in the Galactic plane. As the absolute 
error in the vertical action depends on the sum of the radial 
and vertical actions a large radial action can lead to a large 
relative error in the vertical action. The systematic error 



Table 3. RMS deviations between actions and angles for stars 
in the Geneva- Copenhagen Survey sample estimated with the 
Stackel fitting method and the data from lMcMillanl (|2011rJ ). The 
second column gives the RMS of the expected systematic errors 
from the Stackel fitting method. 



RMS 
difference 



Expected RMS 
error 



AJ i? /10- 4 kpc 2 Myr- 1 


5.4 


4.9 


AJ 2 /10- 4 kpc 2 Myr- 1 


7.6 


4.9 


A^/O.Olrad 


2.6 


0.9 


A<9 /O.Olrad 


1.2 


0.6 


A^/O.Olrad 


7.7 


40.1 



in Z depends on this large relative error. For all the other 
variables the discrepancy between our data and McMillan's 
data seems to be in agreement with the expected systematic 
errors of our method. 



8 CONCLUSIONS 

We have detailed a method for estimating the angle-action 
variables in a general axisymmetric potential given a 6D 
phase- space point. The method is based on locally fitting 
a Stackel potential to the region of the potential the orbit 
probes and then taking advantage of the ease with which we 
may calculate the actions and angles in a Stackel potential. 
We have investigated the systematic errors by producing 
phase-space points of known actions and angles using the 
torus machine and then assessing how well the method can 
reproduce these variables. For a single torus the errors in the 
angles are largest for phase-space points near apsis and the 
errors in the actions are constant (of order a few percent). 
For a collection of tori, chosen to be representative of both 
disc and halo- type tori, the absolute error in the actions is 
found to scale with the sum of the vertical and radial ac- 
tions. The errors in the angles scale with the relative error 
in their corresponding action. We compared the method to 
other methods for estimating the actions in an axisymmet- 
ric potential. The method gives results approximately two 
orders of magnitude more accurate than assuming the poten- 
tial is spherical and performs approximately a single order 
of magnitude better than the adiabatic approximation. 

We have demonstrated that the procedure is suitable 
for most disc and halo-type orbits. The procedure will 
not work for resonant orbits or chaotic orbits. However, 
the occurrence of these orbits in realistic galaxy axisym- 
metric potentials is rare and the great majority of stars 
are on quasi-periodic no n-resonant orbits (|011ongren| [l962; 
iMartinet fe Mavedfl975h . 

We demonstrated the use of the method by application 
to the Geneva- Copenhagen Survey (GCS). As this survey is 
only local, the angle-action space does not reveal much more 
information than velocity space. However, we present angle- 
action coordinates for the peaks of the clumps and streams 
present in the survey and use them to study the relative im- 
pact on estimated angles and actions of observational errors 
and the known systematic errors of the method. We show 
that the observational errors are dominant. 

It is hoped that this method will lead to more 
widespread use of angle-action variables when analysing 
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Table 2. Velocities with respect to the local standard of rest and angle- act ion coordinates for density peaks of known structures in the 
Geneva- Copenhagen Survey. The errors in U and V are 1.7kms~ 1 and 1.6kms — 1 respectively. The errors in the angles and actions are 
presented as a b c where a is the total error and b and c are the contributions from the systematic errors of the Stackel fitting method and 
the errors in the space- velocities respectively, such that a 2 = b 2 ±c 2 . In all cases the observational error dominates the systematic error. 

[//kms- 1 y/kms- 1 J^/kp^Myr- 1 L 2 / kp^Myr" 1 R 



Hercules 


-22.1 


-40.1 


(4.0±0.3g-§ 06 ) x 10" 


2 


1.69 ±0.01 


3.63 ± 0.040;g§| 


(9.6 ± 0.7°;°|) xlO- 2 


Hyades 


-28.2 


-10.8 


(1.2±0.lO : 5 01 ) x 10- 


2 


1.94 ±0.01 


4.37±0.07° : °° 5 


0.118 ±0.007g : oooi 


Pleiades 


-9.8 


-16.3 


(7.3±0.1° : ° 05 ) x 10 - 


3 


1.89 ±0.01 


3.61±0.080 ; ^ 


(4.2±0.70-0° 6 ) x lO- 2 


Coma Berenice 


-2.2 


0.3 


(1.7±1.4° : ° 002 ) x 10 


-4 


2.03 ±0.01 


5.0±l.l?-00° 9 


(9.3±6.4H 008 ) x 10- 3 


Sirius 


14.7 


5.7 


(3.5±0.80 ; o 01 ) x 10- 


3 


2.07 ±0.01 


1.01±0.13°-°° 3 


6.22±O.Ol'o-000° 2 



data, and we intend to release the source code soo iS Whilst 
the GCS can be easily analysed in velocity space, angle- 
action variables should enable us to reveal structures, which 
are more dispersed in phase-space, in larger surveys. 

We have limited the discussion in this paper to ax- 
isymmetric potentials. Whilst for our own spiral galaxy 
the axisymmetric approach may suffice, for analysing el- 
liptical galaxies a triaxial approach must be d eveloped. 
Ther e are also triaxial Stackel potentials (see Ide Zeeuwl 
1985) so it should be possible to expand the approach out- 
lined in this paper to triaxial potentials. The extension 
of the angle-action estimation is simple, but the fitting 
procedure is more compl e x wh en the potential is triaxial. 
Ide Zeeuw fc Lvnden-Belj {l985) discuss how a general tri- 
axial potential may be fitted both locally and globally by a 
Stackel potential. The method for global fitting is the three- 
dimensional generalisation of the method used in this pa- 
per so involves multiple multi-dimensional integrals. Also 
the best choice of coordinate system involves minimising 
the least-square difference with respect to two coordinate 
parameters so a more computationally expensive procedure 
than the simple method used in this paper may be required 
for finding the best coordinate system. 
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APPENDIX A: COMPUTING THE 
ANGLE- ACTION VARIABLES 

The approach taken here as well as the majority of the for- 
mulae have been taken from lEvrel (|2010h . Following on from 
equation J5]) the action J T is given by an integral over a full 
oscillation in r. A full oscillation in A involves integrating 
twice over the interval (Ao,Ai). Ao and Ai are the roots of 
px which may be found by Brent's method using equation 
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(0. There is a complication when calculating J v due to the 
definition of v. v is only uniquely defined for z ^ such that 
a full oscillation in v corresponds to half an oscillation in z. 
Therefore we calculate J v by integrating four times over the 
interval (zy , ^1), where vq = c 2 , as all orbits cross the z = 
plane, and z/i is the root of p u found by Brent's method. The 
actions are given explicitly by 



Jx 



i r Xi 2 r 1 

- / p\d\, J„ = - / p^dzA 



(Al) 



In order to calculate the angle coordinates we use the 
generating function, S, defined as 



S(T,<f>,E,L z ,I 3 ) = S4,+ 



I L z d(j) + ^2 [ Pr'dr' 



(A2) 



This generating function defines the canonical transforma- 
tion between the canonical coordinates (r,<p,p T , L z ) and 
(J T , L Z ,0 T , 0(f)). The angles are now computed as 

_ dS_ _ dS_dE_ dS^dL^ dS_dI^ 

T ~ dJ T ~ dE dJ T + dL z dJ T + dl 3 <9J T ' ^ ' 

for r = A, v. The derivatives of the classical integrals with 
respect to the actions may be found by inverting the 3-by-3 
matrix of the derivatives of the actions with respect to the 
classical integrals. These derivatives are simpler to calculate 
as they follow from equation (|A1[) and the definition of p T 
from equation (|4j): 



dJx = 1 f Xl 
dE 4n J Xn 



dX 



A (A-a 2 )p A : 



dJx 
dL z 



47T 



rX 1 



dX 



(A- 



z px 



dh ^ An 



dA 



/ Ao (X-c?)(X-c^p x 
dJ„ 1 r 1 dv 



J 



8E 2nJ„ a {v-a*W 



dJ„ _ L z 



f" 1 Av 
9/ 3 " 2ttL (v-a 2 )(v- 



(A4) 
(A5) 
(A6) 
(A7) 
(A8) 
(A9) 



The derivatives of the generating function with respect 
to the classical integrals may be calculated in the same spirit 
as 



ds_ _ ^ i r oV 



)Pr 



(A10) 



OS 
dL z 



'■£t{™ (An) 



We note that equation (|A11|) is simply the angle conjugate 
to L z , 0(f). 

With the scheme given above there is a degeneracy in 
U between points in the orbit at ±z. This is simply resolved 
by adding 2ty to U if z < 0. We then must divide U by two 
to ensure 6 U is confined to the interval (0, 2tt). 

As p 2 (r) vanishes at the endpoints of many of these in- 
tegrals, we want to avoid evaluating the integrands at the 
endpoints. We do this by performing a change of variables 
and estimating the integral using a Gauss-Legendre quadra- 
ture scheme. Here we will outline the procedure for calcu- 
lating J\ but the same principle follows for the rest of the 
integrals. We perform a change of variables to 



A = Asintf + A; A = -(A + Ai); A 



1 



(Ai-Ao), (A13) 



' 2 ' 2 > 



(A14) 



such that the integral is now over $ = (- 

i r /2 - 

Jx = - Acostf;p(A(tf))dtf. 

™ J-tt/2 

This integral can now be computed numerically using a 10- 
point Gaussian-Legendre quadrature scheme. 



APPENDIX B: DERIVATION OF BEST-FIT 
STACKEL POTENTIAL FUNCTIONS 

To find the best-fit Stackel potential we must minimise equa- 
tion (fTT]) with respect to the function /. It is useful to con- 
sider minimisation with respect to the two parts of the func- 
tion, /(A) and f(y). This yields 



s: 



dAA(A)[ x (A^)-/(A) + /H]=0, 



dv N(v) [x(A, v) - /(A) + /(i/)] = 0. 



(Bl) 



Rearranging each of these, and noting that A(A) and N(v) 
are normalized over the integration range, we find 



/(A)=x(A)+/ &vN(y)f(y), 



/M = -xH + ^ A+ dAA(A)/(A), 



(B2) 



where the definition of %(r) is given in equation fT4|) . Sub- 
stitution of the expression for f(y) into the expression for 
/(A) we find 



/(A) = x(A) - x 



+ ^ A+ dAA(A)/(A), (B3) 



where x is defined in equation (fT4)) . This equation along with 
IB2I implies that 



L 



i: 



dA A(A)/(A) — / di/JV(i/)/(i/)=£. (B4) 



As we are only constraining the difference [/(A) — f(y)\ we 
are free to choose the values of these integrals as long as 
their difference equals x- We opt for the symmetric choice 

/(A) = x(A) - \x, fW) = -xiy) + \x- (B5) 
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